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We consider a new class of non-Hermitian random matrices, namely the ones which have the form 
of sums of freely independent terms involving unitary matrices. To deal with them, we exploit the 
recently developed quaternion technique. After having derived some general identities describing 
additive properties of unitary matrices, we solve three particular models: CUE plus CUE, CUE 
\Q , plus . . .plus CUE (i. e. the sum of an arbitrary number of CUE matrices), and CUE plus CUE. 

By solution of a given model we mean here to calculate the borderline of the eigenvalues' two- 
dimensional domain, as well as the eigenvalues' density function inside the domain. We confirm 
£SJ , numerically all the results, obtaining very good agreement. 
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I. INTRODUCTION 



Preface 



o 

oo 

O . 

Recently there is growing interest in non-Hermitian random matrix models. One finds them twofold interesting. 
First, non-Hermitian random matrix models are fascinating from the mathematical point of view, because their 
spectra cover two-dimensional and often multiple-connected supports on the whole complex plane, on the contrary to 
Q-r Hermitian ensembles, whose real eigenvalues form cuts on the real axis. This is enough to make most of the methods of 
Hermitian random matrix theory fail. Second, they are ubiquitous in different fields of physics and interdisciplinary 
sciences. It is enough to mention open chaotic scattering 0], spectral properties of Euclidian Dirac operators in 

S the presence of chemical potential [3, the CP-violating angle 9 in QCD non-Hermitian generalizations of the 
Anderson localization in mesoscopic systems Q, modeling of chemical transitions in dissipative systems 5], matrix 
generalizations of multiplicative diffusion processes @ or evolution of spectral curves of non-Hermitian ensembles in 
the context of the growth problem . 

Recently, a new technique of dealing with non-Hermitian random matrix models in the large-./V limit has been 
established ; §| , under the name of quaternion approach, based on free random variables theory [3 . It allows to solve 
non-Hermitian models that acquire the form of sums of "independent" (free, see below) random matrices, in a simple 
algebraic way. 

In this paper we shall concentrate on a new class of non-Hermitian sums of random matrices and apply the 
quaternion technique to solve them. These will be models having unitary matrices as summands. Unitary random 
matrices find applications e. g. in chaotic scattering [Toj . conductance in mesoscopic systems [TTI | or periodically 
driven quantum systems |12| . The most important ensembles of this type are the Dyson models [l3| . i. e. the 
so-called COE, CUE and CSE; see also for a detailed discussion. We shall pick CUE as our basic example. For 
a variety of issues concerning unitary models we refer to • 

In other words, we shall try to add unitary random matrices. This may seem quite odd because what one usually 
does with unitary matrices is to multiply them since the product of unitary matrices remains unitary, which is not 
the case when adding them. Moreover, any possible applications of such sums of unitary random matrices are still 
unclear. Nevertheless, wc think that the problem is interesting at least from the mathematical point of view and its 
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solution reveals how easy it is to treat non-Hermitian sums of random matrices using the quaternion technique, even 
if these sums appear very involved. 

The rest of sec. [I] is devoted to short summary of the quaternion technique. In sec. ^ we apply the quaternion 
method to the case of sums of free unitary random matrices, deriving a number of useful formulae. Sec. IIIII llVl and Ivl 
contain solutions of three examples: the sum of two free CUE matrices, the sum of an arbitrary (and eventually also 
infinite) number of free CUE matrices and the sum of the free CUE and CUE matrices. All the results are confirmed 
by numerical computations. Finally, sec . I VII summarizes the paper and points to some future possibilities of exploring 
the subject. 



B. Summary of Quaternion Approach 



Let us here very briefly recall the foundations of the quaternion approach; details and further references can be 
found in 

Free random variables. The theoretical background of the quaternion approach contains the free random variables 
calculus (FRV) |9( , which is a generalization of classical probability calculus to the case where random variables 
become non-commutative. In such a non-commutative setting it is possible to define a pertinent extension of 
the notion of independence, called freeness, which allows to lift the ideas known from classical probability and 
built upon the concept of independence onto the level of non-commutative random variables; these ideas are e. 
g. additive and multiplicative convolution, central limit theorems, etc. 

We shall deal with one of them, the problem of adding two free non-commutative random variables, a\ and a2- 
On the classical level this is readily solved by introducing the characteristic function g a (z) = 53 n >0 ^7rr z " = 
(e 2a ), which collects all the moments in a single generating function; and its logarithm satisfies the addition 
law, In g ai+a2 (z) = \ng ai (z) + lng a2 (z), for two independent commutative random variables a\ and 02- Hence 
in classical probability it is just logarithm that makes an additive function from the generating function of the 
moments. On the non-commutative level the task is much more involved but the final algorithm is analogous: 
The moments are gathered in another generating function, called the Green's function, 

n>0 N ' 

and what needs to be done in order to obtain an additive function under the addition of two free non-commutative 
random variables is not to take logarithm but to invert it functionally, 

Ga(B a (z)) = B a {G a {z)) = z, (2) 

where this new object is called the Blue's function and obeys the addition law, 

B ai+a2 {z) = B ai (z) + B a2 (z)--. (3) 

z 

Hermitian random matrix theory. Random matrix theory (RMT) in the large-TV limit, which deals with large 
N x N random matrices, is a particular instance of the FRV calculus, with the following expectation value map, 

(X)r.n, = (^ Trf ) i > ( 4 ) 

where (. . .) c ], is the classical expectation value. We shall focus on the eigenvalues of a given random matrix 
model. If we restrict to Hermitian matrices, call H one of them, the eigenvalues Ai, . . . , Aj\r are real and their 
density is defined through 



PH 



this is the crucial object that we aim to find. It turns out that the same Green's function which is so useful in 
the addition problem in general FRV, and which here acquires the form 

G H (z) = ±(Tr 1 - \ , (6) 
N \ zl N - i27 . 
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(ljv is the N x N unit matrix) serves perfectly also to investigate the eigenvalues' density: it is a meromorphic 
function with poles at the eigenvalues on the real axis, which in the large-TV limit become cuts, and its behaviour 
close to these cuts provides the desired eigenvalues' density, 

p H (X) = -- lim ImG H (A + ie). (7) 

Non Hermitian random matrix theory. When we release the symmetry requirement that the matrices are Hcr- 
mitian, the theory complicates significantly, which is mainly due to the fact that eigenvalues depart from the 
real line and become complex in general, forming in the large-iV limit two-dimensional domains instead of cuts. 
This deprives the Green's function of its usefulness as it is singular inside the eigenvalues' domains and can no 
longer provide the eigenvalues' density px{z, z) of a non-Hermitian random matrix X. The idea is to regularize 
the Green's function, 

where we put an additional factor of (zl N — X^) both in the numerator and denominator, and moreover the 
denominator is regularized by e 2 lpj with e — > + (we shall skip the limit symbol henceforth). This is called 
the non-holomorphic Green's function. It is defined everywhere, also inside the eigenvalues' domains, as the 
regularized denominator is positive definite on the whole complex plane; however, it is non-holomorphic. Not 
only have we now gained access to eigenvalues' domains but we can also reproduce the eigenvalues' density from 
the non-holomorphic Green's function via the simple formula 

px(z,z) = -d,G x (z,z). (9) 

Even though the non-holomorphic Green's function is so meaningful, it seems to be very involved to evaluate 
due to the quadratic structure in its denominator instead of the linear one for the complex Green's function. 
This obstacle can be overcome by introducing yet a new object, the matrix-valued Green's function, which is a 
2x2 matrix defined as 



2N x 2N 




where the block-trace operation is 



and for short 



. . A B \ ( TrA TrB , 

^ D) 2Nx2N ={TrCTrD) ("> 



Ze = ( I 'I ) - (12) 



ie z 



2x2 



•v" ( A A . ) • (I*) 



2JVX2JV 



This matrix looks as the usual complex Green's function but lifted to the level of 2 x 2 matrices and this 
similarity implies that some methods of Hermitian RMT used to derive the complex Green's function can be 
applied also in non-Hermitian RMT for the matrix-valued Green's function. Moreover, its 11-element equals 
precisely the non-holomorphic Green's function, 

G x (z,z) = [g x (z,z)] lu (14) 

which in this way becomes accessible. 

The rest of the matrix-valued Green's function, i. e. its off-diagonal elements, may seem to carry no information, 
but actually they are important as well: their product 



C x (z,z) = [Gx(z, z)} 12 [gx(z,z)] 21 , 



(15) 
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evaluates the correlator between left and right eigenvectors of X, 

ImiLiXRil^S^iz-Xi)) =--C x {z,z). (16) 



N 

£( 

i=l 



TV 



We shall exploit one of the properties of this object, namely the fact that it vanishes on the borderline of the 
eigenvalues' domains and remains zero outside: hence a very important geometric characteristics of a given 
non-Hcrmitian model, the equation of the borderline, is readily obtained from the correlator, 

C x (z,z) = 0. (17) 

Therefore we shall be interested in the whole matrix-valued Green's function as it provides us with both the 
eigenvalues' density and the shape of the eigenvalues' domains. 

Quaternion formalism. This new approach makes use of insights from both the FRV calculus and non-Hcrmitian 
RMT. The matrix-valued Green's function is generalized to a quaternion function of a quaternion variable, 
where we recall that a quaternion is a 2 x 2 matrix defined by two complex numbers a and b via 

by replacing Z e with an arbitrary quaternion Q, which gives the quaternion Green's function, 

g x (Q) = l.( h Tr —-l vn ) . (19) 

cl 



N\ Q<E>1 N -X D 



For Q = Z e we regain the former meaning together with its usefulness, but the complete functional dependence 
on Q proves to be crucial. Indeed, it allows, following the FRV analogy, to invert Gx(Q) functionally to get the 
quaternion Blue 's function, 

Gx(B x (Q)) = B x {Qx(Q)) = Q. (20) 

This new object might have turned out to be spurious as it generalizes just some notion from non-Hermitian 
RMT, whereas the Blue's function is a FRV concept. But it may be shown that in a general setting of the 
so-called FRV calculus with amalgamation, which is an extension of usual FRV, the quaternion Blue's function 
plays a fully analogous role to the usual Blue's function, namely it obeys the quaternion addition law, 

B Xl+ x 2 (Q)=B Xl (Q) + Bx 2 {Q)-^ (21) 

for two free random matrices X% and X^ , Hermitian or non-Hermitian. Therefore the addition problem in non- 
Hcrmitian RMT is solved: it relies on the knowledge of the quaternion Green's functions of both summands, X\ 
and Xi, inverting them functionally, exploiting the quaternion addition law to get the quaternion Blue's function 
of X\+Xi, and finally inverting it functionally at the point of Q = diag(z, z) to obtain the matrix-valued Green's 
function of X\ + X2 , i. e. solving 

B Xl +x 2 {Qx x +x 2 (z, z)) = diag(z, z), (22) 

where we can set e = as the regularization procedure happens to be entirely encoded in the functional inversion. 
This technique is called the quaternion addition algorithm. 



II. QUATERNION FORMALISM FOR UNITARY RANDOM MATRICES 



A. Introduction 



Let us start with the following observation, which continues the considerations from the basic paper gj. We have 
seen that it is the quaternion Green's function Qx{Q) EH which is the basic object to deal with a non-Hermitian 
random matrix model X within the quaternion formalism. For instance, to apply the quaternion addition algorithm 
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we need to know the quaternion Green's functions of both summands. Hence it becomes crucial to have a means 
of computing this object for some classes of random matrices. Now it is observed that there are cases where the 
quaternion Green's function can be quite easily calculated explicitly. More precisely, there are cases where it can be 
reduced to the complex Green's function Gx{z) ©, 

G x (z) — g x (Q); (23) 

and this passage is called the Hermitization procedure. In other words, there can be identified instances of X such that 
the whole information contained in the quaternion Green's function follows from a much simpler object, the complex 
Green's function. 

To see how such an idea could arise, let us note that Gx(Q) is the averaged block trace of the inversion of a 2N x 2N 
matrix, and let us simply perform this matrix inversion, 

Gx{Q) = - (bTr (Q u - X^Y 1 ) = i /bTr ( c1n ~ X id " ljv vt V * \ = 
{ ' N\ y ^ > / N\ \ idl N cljv-Xt J 2Nx2N / 




cX + (\c\ 

— id cl/ 



XXi -cXI -cX+(\c\ 2 + \d\ 2 )l N XXi -cXI -cX+(^c\ 2 + \d\ 2 )l N J \ ^4) 



XtX-cXI -cX + (\c\ 2 + \d\ 2 )l N Xi X-cX1 -cX + (\c\ 2 + \d\ 2 )l N 

where now we denote the quaternion Q by 



2Nx2Nl 



2x2 



to agree with our later conventions. 

Now every block of the above 2N x 2N matrix is seen to be a rational function of X and X'. The denominators 
here are quadratic, which is in contrast with the linear denominator featuring in the complex Green's function JJjJ, 
and this mainly hinders efficient computation of Qx{Q)- However, if there is a symmetry restriction on the matrix X 
which assumes the following form, 

X T = rational function of X, (26) 

then every block of the matrix in question becomes actually a rational function only of X. If so, it can be expanded 
in simple fractions, i. e. functions of the form l/(sljv — X), for some complex numbers s; this is nothing but the 
structure appearing in Gx(s), hence the full quaternion Green's function seems likely to be expressed solely in terms 
of the complex Green's function, thus proving the Hermitization procedure (|23|l for X . 

The symmetry constraint ()26|l . even though at first it looks quite fanciful, is not of minor importance. There are at 
least its two important realizations, 

i? f = H (Hcrmitian), C/ 1 = U^ 1 (unitary). (27) 

In both these cases, an appropriate Hermitization procedure should be applicable, thus reducing a complicated 
quaternion object to a much simpler comlpex one. 

In the Hermitian case, X — H, the quaternion Green's function has been Hermitized in the basic paper and let 
us just quote the resulting formula for further comparison, 

QH{Q)=lH{q,q)l2-i H {d,Q)Q ] , (28) 
where and j' H are two scalar functions depending only on the eigenvalues q, q of Q, and given by 

, _ qG H (q) - qG H (q) ,, _ G H {q) - G H (q) 
lH{q, q) = = , 7h(?> q) = = • (29) 

q-q q-q 

This allows to have the full quaternion object Qh{Q) once the complex one, Gh(z), is known. 
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B. Quaternion Green's Function for Unitary Random Matrix 



The same idea can be now applied also in the unitary case, X — U. The constraint U> — U^ 1 shall be substituted 
to (|24|l . and the fractions expanded in simple fractions. Due to the quaternion structure of the result, it is sufficient 
to write only e. g. the left-upper and left-lower elements, the other (depicted by stars) simply follow from these two 
via appropriate complex conjugation, compare (|25J) . 

i / / ci N -y- L -id 
Qu{Q) = ^ (bTr[ ^- cl Tu~ 5U 



gl N -cU 1 —cll gl N -cU 1 -cU / 2Nx2N 




2N x 2N ' 



iu{Q)-\iu{Q) ★ 



(30) 



where we introduce the following notation, 



lr .,_u 1 G u (u 1 )-u 2 G u {u 2 ) Guiu^ - Gu(u 2 ) , , 

lu{Q) = , Ju\Q) = ) 

m -u 2 Ul- u 2 

where u\ t 2 are the solutions of the quadratic equation cu 2 — gu + c = 0, 

«i,a = ^(ff± V^_4|c|2), (32) 

where for short 

.g^| C | 2 + |d| 2 + l. (33) 

In this simple way we are given the quaternion Green's function for any unitary random matrix model U , expressed 
via its complex Green's function. 

Note similarities and differences between (|30|l . I|31|l in the unitary case and l|28|) . (|29|l in the Hcrmitian one. In 
particular, in the Hermitian case the gamma functions depend on the eigenvalues q, q of Q, so they are invariant under 
similarity transformations of Q, whereas in the unitary case we have a more involved and non-rotationally invariant 
dependence of the gamma coefficients on u\^ 2 . 

Note that u\u 2 = mua = 1, hence we can denote 

u = ui, and consequently u 2 = — , (34) 



which gives, due to the formula Gu{z) = i (l — -|G(y(i)), 



uGu^+uGuju)-! Gu^+u^uGuju)-!) 

iu(Q) = 1 , iu(Q) = 1 • ( 35 ) 



U — T U — ■= 



Thus we have expressed the gamma functions through the single complex variable u instead of two mutually dependent 
variables u\ }2 . 

Let us note that we have obtained the basic expression l|30|l under some slight assumptions, 

• ui u 2 , i. e. \c\ 1 or d ^ 0. 
From now on we have to consider these two special cases separately from the generic case. For c — we get 

1 ( mUJ id \ r „ fi \ 
\d\* + l\ id mu,i) 2x2 > W 

where mu,i — -^-(TrC/) is the first moment of U. For |c| = 1 and d — immediately 




^(diag(c, c)) = diag(G c/ (c), Gu(c)); (37) 

we have to pay attention when using this formula since here |c| = 1 and unitary random matrices have their eigenvalues 
exactly on the unit circle, so that Gjj{c) is given by a divergent integral, which however makes sense as a certain limit. 
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C. Example: Quaternion Green's Function for CUE 

Let us now consider a particular instance of a unitary random matrix, namely the circular unitary ensemble (CUE) , 
which means 

P(U) = const, (38) 

or in other words, the eigenvalues' density reads 

Pu(z 7 z) — — on the unit circle. (39) 

The complex Green's function is derived easily, 

1 f 271 d9 i f 2 " de~ 10 i f dt f Res^i^, for \z\ > 1 



i f de i r 

2?r J z- e tV 2tt J 



2tt J ze' 18 - 1 2-7T J C (o,i) zt - 1 [0, for \z\ < 1 



i, for \z\ > 1 
0, for \z\ < 1 



(40) 



therefore 



In order to calculate the quaternion Green's function, we have to find the gamma functions. Since there is always 
\u\ > 1 due to \u\ > 1 <^> (|c| - l) 2 + \d\ 2 > 0, we have 

C C QC 

lu{Q) = — r = = jz , iu{Q) = -rr + ~ — / 9 

Vd 2 ~ 4|c| 2 2c 2c^ 3 2 -4|c| 2 

[^(Q)]u = j-+ ' c|2 ,~ M|2 ~, 1 2 , [MQ)]2i = ~ , = , 2 - (41) 

2c 2 Cv / ff 2 -4|c| 2 ^5 2 -4|c| 2 

These constitute the quaternion Green's function for the CUE random matrix. 
For the special case c = 0, since the first moment of the CUE random matrix vanishes, 

*(U?)J-TFM(i?L- 



D. Quaternion Blue's Function for Unitary Random Matrix 

We have seen that it is the quaternion Blue's function Bx(Q) l|2U|) . i. e. the functional inverse of the quaternion 
Green's function, which is the basic object entering the quaternion addition law (|21|l . allowing us to add free random 
matrices. 

Here we aim to invert functionally GuiQ) <|30[) . The resulting quaternion Blue's function is a quaternion denoted 
by 

*«>-(«?),„■ <«> 

where we return to the notation from sec. [I] for the quaternion Q, 

(44) 



a ib 
ib a 



2x2 



It is provided by solving the two equations 



lu (B u (Q))-^j' u (Bu(Q)) = a, ~lu(Bu(Q))=b, (45) 
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with two complex unknowns, c and d, assuming the knowledge of a and 6, i. e. the complex coefficients of the 
quaternion Q. These equations are valid only for the generic case of c ^ and (|c| ^ 1 or ^ 0). We shall not solve 
them as we shall see that we actaully do not need an explicit solution, barely having equations is enough. 

If c = then the equations are 



which has a solution only while 



which is then of the form 



mu,i = «, = b, 



0<|6|<i, a = mUJ, (46) 



„, ,- ' + V^W (47) 



E. Example: Quaternion Blue's Function for CUE 

In the CUE case we know Gu(z) explicitly, so we can also explicitly write the equations (|45(l . 

' ■ II'- W- 1 - (48) 

This time we are also not going to solve this equations. The reason is that they will soon be used together with the 
quaternion addition law, which imposes some additional constraints on our unknowns c and d, and it will turn out 
that it is easier to solve all these equations together. 

Let us however simplify them a little. Let us assume first the generic case of c ^ 0, a ^ 0. We see from the first of 
these equations that 

ttEacel, (49) 

and we rewrite the second one as 



d = -b/3, /3 = vV -4|c| 2 e M+. (50) 
Hence we are left with the following two equations, 



where 



p=^ 2 - 4 w (51) 

2a- 1=W—LL± , (52) 

P 

g = ^ + H 2 /? 2 + l, (53) 
\<A 

for two real unknowns, a and [3. 

If c 7^ it however may turn out that a — 0, and so the derivation does not hold. Instead, we get after a short 
calculation, 

|c| = l, d = 0, 

which is impossible to be reached from a = due to 1)37(1 . 
If c = 0, one must have a = and < |6| < |, and then d is given by l|47|l . 
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III. CUE PLUS CUE MODEL 

We have established foundations of the quaternion approach in the case of unitary random matrices. Now let us 
move to a number of random matrix models which have the form of a sum of matrcies such as CUE or GUE. 



A. Introduction 

As the first example we consider the model 

X = U 1 + U 2 , (54) 



where L/1,2 are two free CUE random matrices. 
The quaternion addition law l(21|) reads 



i. e. denoting 



B Ul+ u 2 (Q) = B Vl (Q) + B U2 (Q) - i = 2Bu(Q) - 1 



B Ul+Ua (Q)=( w'c ) ■ ( 55 ) 



and exploiting our previous notation, 



C = 2c- . ., ° D = 2d+ b (56) 

|a| 2 + |6| 2 |a| 2 + |6| 2 



To solve the model (i. e. to find the eigenvalues' density and the borderline of the eigenvalues' domain) we need to 
consider the matrix-valued Green's function l|lUfl . i- e. to solve the equation (|22H . 

Bu!+u 3 (Q) = diag(z,z), 
with respect to Q, where z is a given complex number. This reads 

z = 2c ~WTW' (57) 
= 2d+ WTW (58) 

Our aim is now to solve these equations with respect to a and b, where c and d are given by 1149(1 . (|50|l . I|51() . I|52|l . 
1)53(1. Let us note that here a^O, because otherwise c = which would contradict 1 (57(1 . 



B. Solution 

Let us write down all the equations we have to solve. First, 1(57(1 multiplied by a shows that 

7 = za S R, (59) 

which in particular means that 

M 2 = ]i. (60) 
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so that (|57|l assumes the form 

" = 2 ' 2( 7 2 + |z| 2 |6| 2 ) 



9 

7 , 7 



" = o + o/-2^U|2IM2V ( 61 ) 



Second, together with (|5T]l acquires the form 



^ 2( 7 2 + |l| 2 |6| 2 )- (62) 



Third, let us recall j52)l. (f53jl in the new setup, 



«. e. 



/?=J.9 2 -4^, (63) 



^-|6| 2 /3 2 -l 

2a-l = ^ Ui , (64) 

P 

+ + (65) 

7 

To sum up, we aim to solve five real equations, (|5T)) . (JHU, @, l|t)5fl. for five real unknowns, 7, \b\ 2 , a, /3, g, 

among which we are specifically interested in 7 and \b\ 2 . 

To exectute this plan, we substitute (|61|l and l|62|) into the remaining equations. First, into Ij64(l . which simplifies 
surprisingly into 

2 _ l Z l I |2|i|2 



Second, to 

|z| 2 



.9 = 2 + 2 7 ( 1 
Third, to (|63[) . which gives a half of the solution, 



„ _2(2-M 2 ) 



I 6 ' = l) ! 2 2 ' ^ 

(4- |z| 2 ) 



so that also the second half, 



4-\z 

This can be rewritten in terms of the basic quantities 



7 = ^-02- (67) 



G Ul+U2 (z,z) = — -, (68) 



4-U 



2(2- \z\ 2 ) 

C Ul+U2 (z,z)= ) \:> . (69) 
(4- \z\ 2 ) 



Therefore, the eigenvalue's density is 



p Ul+ u 2 (z,z) = — —2, (70) 

7T (4- \z\ 2 ) 

which can be easily checked to be normalized to 1, whereas the borderline's equation l|17f) reads 

1*1 = A (71) 
which means that the eigenvalues fill the centered circle of radius y/2. 
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C. Numerical Confirmation 



Let us numerically confirm the above results. We have drawn 100 pairs of 200 x 200 unitary matrices from the 
uniform distribution (CUE), added them, and diagonalized the sum. 

The left figure shows positions of these 20000 eigenvalues on the complex plane as well as the theoretical borderline 
1)71(1. The right figure presents the radial section of the eigenvalues' density plot (for it is circularly symmetric); the 
solid curve stands for the theoretical result l|70|l . and there is a numerical histogram under the curve. 




We are convinced of the perfect agreement between theory and numerics. 

Some slight deviations are due to a finite (N = 200) size of matrices. In the origin there is also a strange peak 
visible, which we suspect to be a numerical artefact. 



IV. CUE PLUS . . . PLUS CUE MODEL 



A. Introduction 



As the second model, let us consider a slight generalization of the CUE plus CUE model, namely the sum of M > 2 
free CUE random matrices f/j, 

X = U 1 + ... + U M . (72) 

In particular, we shall eventually multiply A by a constant depending on M (1/yM, specifically) to ensure that the 
result has a correct behaviour for M — > oo, and we shall investigate this limit. This will be free additive unitary 
diffusion. 

The quaternion addition law l|21|) states that 

Bu 1 +...+u M (Q) = MBu(Q) - ^— 



c = Mc ~" {M ~ l) WTW' D = Md + {M ~ l) WTW- (73) 

Again, the basic equations are 

z = Mc -w-VwTW' (74) 
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Q = Md+(M-l) b (75) 



Again, only the generic case of a ^ needs to be considered. 



B. Solution 



Let us write equations to solve. First, (|?4l multiplied by a gives 

7 ee za £ R, (76) 

with J71J being then 

1 / 1 \ 7 2 

a= M 1+ \ l -M) ^TWW (77) 

Second, (|75|l together with (|50|l assume the form 

={ l -^)^rww- (78) 

Third, the basic equations l|51[) . (|52[) . (|53[) remain identical as in the CUE plus CUE case, 



^J.9 2 -4^, (79) 



|6| 2 /3 2 

2«-l= ^ , (80) 



q 2 |z| 2 

7 2 



+ |6| 2 /3 2 + l. (81) 



We get thus to solve five real equations, lt77|l . f75|) . (|7TTfl . JSUJ), (|B*T)l for five real variables, 7, |6| 2 , a, /3, 5, among 
which we are specifically interested in 7 and \b\ 2 . 

The equations are almost identical to those for the CUE plus CUE model. Hence proceeding analogously, 

z(M - 1) 

Gu 1 +...+u M {z,z) = M2 _ (82) 

M(M - 1) (M - |z| 2 ) 

- k[/ 1+ ...+i/ M (2, zj - 777; —r-2 • K H6 ) 

(M — \z\ A ) 

For M = 2 we regain the results of the previous section. 
The eigenvalues' density thus reads 

, _ N M 2 (M-1) . . 

Pu ^... +UM {z,z) = ^_^ r (84) 

normalized to 1, and the borderline's equation is 

\z\ = VM, (85) 



which means that the eigenvalues fill the centered circle of radius v M. 
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C. Free Additive Unitary Diffusion 

The above solution is valid for any M > 2. Now let us consider the limit 

M -> oo, (86) 

and investigate the behaviour of the eigenvalues' density and borderline's equation when M grows to infinity. 

Obviously, this cannot be done for the model considered above, i. e. U\ + . . . + Um, since it has no proper large-M 
limit; the eigenvalues' density (|84J) tends to zero and the radius of the bordering circle l|85(l grows to infinity. Since 
we have the scaling relations 

I / z z \ I / z z\ ! / z z\ 

G k x(z,z) = -Gx[^,-j:), C kX {z,z) = ^C x {-,-) , p kX (z,z) = ^ Px (87) 



for real fc, we see that we need to choose k — r^/VM, where is an arbitrary real and positive constant, which 
gives, still for finite M, 

G u 1 +...+u M (z, z) — — -j M \ r , (88) 

30 -/M r 2 ( 1 _ \ z \ \ 

r (■? -A — ^ ~ ^) ( r °° Z (RQ\ 
— C Ul + ...+u M [z, z) — — , (89) 



rl II 



therefore the eigenvalues' density reads 



p r _u 1 ±_^u M _(z,z) = , (90) 



and the borderline is described by 



(91) 



Now we are ready to make the M — > oo limit. The non-holomorphic Green's function and the eigenvectors' 
correlator tend to 

lim G u 1+ ...+u M (z,z) = - r , (92) 

- lim C a 1± ^ L (z,z)= r ^~ |z|2 , (93) 

hence the limiting eigenvalues' density reads 

lim p u 1 +...+u M (z, z) — — — — const, (94) 

whereas the borderline is independent of M and remains the centered circle of radius r M . In other words, the 
eigenvalues are uniformly distributed inside the centered circle of radius r^. Note that this distribution is the same 
as the eigenvalues' distribution of the Girko-Ginibre model 0|. We can regard the described property as a certain 
kind of central limit theorem. 
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D. Numerical Confirmation 



We shall confirm the results for a few values of M, namely for M = 3, 5, 10. For each of them we have drawn 100 
times M unitary matrices of size 200 x 200 from the uniform distribution (CUE), added them, normalized through 
dividing by vM, and diagonalized the result. We have chosen = 1. 

The left column of the figures shows, for each of the above values of M, positions of these 20000 eigenvalues 
on the complex plane as well as the theoretical borderline 1|85|) ■ The right column presents the radial sections of 
the eigenvalues' density plots (again, they are circularly symmetric); the solid curves are theoretical (jHU, and the 
histograms are numerical. 




numeric^ eigenv^ues 



0.5 
0.45 



numerics - 
theory - 
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Again, the complete agreement between theory and numerics is found. 

We also see how the density converges into the uniform distribution on the unit circle, in agreement with l|90|) and 

V. CUE PLUS GUE MODEL 
A. Introduction 



Let us now solve the following model, 



X = U+pH, 



(95) 



where U is the CUE random matrix, H is the GUE random matrix, and they are free, whereas p is a given real 
constant, which may be assumed to be positive, for we shall see that p appears everywhere only as p 2 . 

The quaternion addition law l|21|) is now 

B u+v h{Q) = Bu(Q) + B p h(Q) - ^ = Bu(Q) + p 2 Q, 

where we have used the formula for the quaternion Blue's function for GUE and the scaling law, B p h{Q) = pBh(pQ) = 
p 2 Q + 1/Q, i. e. 



The basic equations are therefore 



C = c + p 2 a, D = d + p 2 b. 



z = c + p a, 



(96) 



(97) 



= d + p 2 b. (98) 
Once again, we are only to consider the generic case of a ^ 0, because otherwise c — 0, which would contradict (|97|) . 



B. Solution 



First, together with g^J show that 

a = za — p 2 a 2 . 



(99) 



16 



Therefore we break a into its real and imaginary parts, 

a = oj + iuj', (100) 

which changes Ij99(l into two real conditions, 



a = xuj - yui' - p 2 uj 2 + p 2 u}' 2 , 



= xuj' + yuj — 2p 2 ujuj' . 
The second one can be used to express uj' through u>, 



yuj 



2p 2 uj — X ' 

and hence also a through uj, 



(101) 



« = a*, - pW + p 2 (-^-) 2 - -p-. (102) 
\ 2p z uj — x ) 2p z uj — x 



If we compute uj, this will give also a. 
Let us note that we can express \a\ 2 through uj, 



MWWW + fj-^-) 2 . ,io 3 ) 



One can also check the simple identity, 



^>=--P 2 > (104) 
or u) 



which will soon become useful. 
Moreover, l|98|l states together with (|50|l that (3 is just a constant, 



d = v\ i- e. ^g 2 -A^=p 2 . (105) 

Let us first inspect consequences of these simplifications for l|52|) . which expresses, after some manipulations, \b\ 2 
through uj and a, 



\ b \p-{~-Zp) a +P- 1 = 



or only through uj if we exploit l|102|) 

1 



(x - 2pujf 



(x 2 (p 2 - 1 + x 2 + y 2 ) + jkc(4(1 - p 2 ) - (5 + ip)x 2 - (1 + 3p)y 2 )uj+ 



+p 2 {4{p 2 - 1) + (8 + 15p)x 2 + 3py 2 )uj 2 - 4p 3 (l + 6p)xuj 3 + 12p'V) . (106) 
If we compute uj, this will give also (3. From l|53|) we get 



x 



g = 2[--2p 2 )a+p 2 , (107) 



and so the combination 

5 2 -4— -/3 2 = 4 (--2p 2 ) a 2 +A{-(p 2 -l)+p 2 {l-2p 
or \uj / \uj 



a. 
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which gives due to i|51[l an equation for u>, 

(1 - 2 P A 2 a = -(1-p 2 ) + p 2 (2p 2 - 1), (108) 

\ to / CO 

i. e., after using the explicit form of a l|102|) . 

{p 2 - l)x + x 3 + xy 2 +p 2 (l - 2p 2 - 5x 2 - y 2 )uj + 8p 4 xuj 2 - 4p 6 u; 3 = 0, (109) 

which is a desired third order (Cardano-type) equation for lu. This is the solution, because we have found the equation 
satisfied by ui and expressed a and |6| 2 through u>. 

Let us however remind that all the calculations are done under the assumption of the generic case of c ^ and 
a 0. Let us thus investigate whether there exists the limit of 

a — > 0, which corresponds to ui — > (110) 

of our solutions. Immediately we see from (|106|l that \b\ 2 p A tends to 

|&|V ->p 2 -1 + x 2 +y 2 , (111) 

hence the equation l|109fl for ui simplifies into 

P 2 - 1 + x 2 + y 2 = 0, 

which means that there exists a solution of our basic set of equations which cannot be approached from the generic 
case, and which reads 

6 = 0, (112) 

and 

x 2 +y 2 = l-p 2 . (113) 

This solution exists only when (|113f) has solutions for x and y, i. e. for 

p<l. (114) 

Let us interpret this solution. The condition (|112fl of vanishing of b means <|17[) that this solution is valid on the 
borderline of the eigenvalues' domain. Therefore (|113f) says that for p < 1 at least a part of the eigenvalues' domain 
is a circle of radius \/\ — p 2 . We will see that this is indeed a part of the borderline, namely the internal boundary, 
whereas the external boundary is an ellipse. We may regard this limiting procedure as investigating the so-called 
holomorphic limit |8|], since we approach the borderline of the eigenvalues' domain. 

To summarize, the solution is given by the Cardano-type equation (|109fl for id, which provides a via (|10Uf) and 
(UnU, as well as |6| 2 via l(TTOI) . 

&u+ p h{x, y) + i 2 _ , (115) 

-C u+pH (x, y) = p ^ x ^ 2pu)2 V - 1 + x 2 + V 2 ) +Ps(4(l - P 2 ) - (5 + 3p)z 2 - (1 + 3p)y 2 )io+ 

+p 2 {4:{p 2 - 1) + (8 + I5p)x 2 + 3py 2 )u 2 - 4p 3 (l + 6 P )xlu 3 + 12p 5 w 4 ) . (116) 
The eigenvalues' density follows thus immediately, 

Pu+ p h{x, y) = (d x + id y ) (ui + iu)') = (d x u; - d y u;') = . .. 
Ztt Ztt 

where one can check that d x u>' + d y oj = 0, so that the imaginary part vanishes, as expected. Exploiting l|102|l we have 
further 
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Differentiating now (|109fl with respect to x and y we get 

1 - p 2 - 3x 2 ~ y 2 + 10 P 2 xuj - 8p 4 w 2 



d x LO 



dyUi 



p 2 (l — 2p 2 — bx 2 — y 2 ) + lQp A xuj — \2p 6 Lu 2 ' 
2y(p 2 uj — x) 



y p 2 {l - 2p 2 - 5x 2 - y 2 ) + 16p 4 xuj - 12p 6 u 2 ' 
therefore finally 

pu+ P H(x,y) = 

= (x 2 (l - p 2 ) - 3x 4 - 3x 2 y 2 + p 2 x(-3 + 2p 2 + Ylx 2 + 5y 2 )uj + 2p 4 (l - 17a; 2 - y 2 )u? + 28p £ W - 8p 8 w 4 ) / 
/2tt (-p 2 (a; 2 (2p 2 - 1) + 5x 4 + x 2 y 2 ) - ^p A x{\ - 2p 2 - 9x 2 - y 2 )to~ 

-4p 6 (-l + 2p 2 + 24t 2 + y 2 )uo 2 + 112 P 8 xuj 3 - 48p 10 w 4 ) , (117) 

where u) is given by (|1U9|I . This is the desired eigenvalues' density in a parametric form. 
We need also the borderline's equation, which is \b\ 2 = 0, i. e. from l|106(l . 

1 -p 2 



z-3p 2 ' 

which put into (|108fl gives a surprisingly simple solution 



1 + 2p 2 ' 

this value of to (valid when on the borderline) should be substituted to l|l(J6|l to get the borderline's equation, 

' r * 2 + (l+ P V = l, (118) 



(l + 2p 2 ) 2 ' 

which is an ellipse. If p < 1, this is a part of the borderline, because we know that the borderline can be reached also 
via the a — * limit, which has given the circle 

x 2 +y 2 = l-p 2 . (119) 

To summarize, if p < 1, the eigenvalues of the U +pH model are placed within the domain with two boundaries, the 
external one, which is the ellipse 1)118(1 and the internal one which is the circle l(119|l : for p = 1 the internal circle 
reduces to the point (0,0), so to vanish for p > 1, and the eigenvalues fill the whole ellipse l|118|) . We have therefore 
a phenomenon at p — 1 which may be called a topological phase transition. 



C. Numerical Confirmation 



The results are confirmed for a few values of p, namely for p = 0.5, 0.75, 1, 2. For each of them we have drawn 50 
unitary matrices of size 100 x 100 from the uniform distribution (CUE), as well as 50 Hermitian matrices of the same 
size from the Gaussian distribution (GUE), added them with appropriate p, and diagonalized the sum. 

The first four figures show, for each of the above values of p, positions of these 5000 eigenvalues on the complex 
plane as well as the theoretical borderlines 1(118(1 and (for p < 1) 1(119(1 . We see in particular that the internal circle 
decreases as p increases, and vanishes completely for p — 1 (topological phase transition), so that for p > 1 there 
remains only the external ellipse. 
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-1.5 -1 -0.5 0.5 1 x -2 -1.5 -1 -0.5 0.5 1 1.5 




-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 x -5-4-3-2-101234 



Experimental points lie exactly inside the theoretically predicted domains. 

To make our figures even more expressive when investingating the eigenvalues' density, we have increases the statistic 
of matrices, namely we have drawn 100 sets of 200 x 200 CUEs and GUEs, added them appropriately, and diagonalized 
the sum. We have repeated it for p = 0.5, 0.75, 1, 2. 

The left column of figures below contains numerical histograms in three dimensions of the density, based on 20000 
experimental points each time, whereas on the right we see the theoretical prediction (|117fl . 
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We recognize the plots to be similar. It is striking that such a complicated expression as (|117|l is confirmed 
numerically in such a good way. Slight differences between the detailed shapes of the respective figures are due to 
numerical obstacles in plotting three dimensional numerical histograms, as well as solving numerically the equation 
(|109fl . choosing the proper solution, and inserting it into the expression <|117ll for the eigenvalues' density. (In 
particular, the two holes in the theoretical figure for p = 0.5 are artificial and mean that the plot is relatively thin 
there.) 



VI. SUMMARY AND PROSPECTS 



In this paper we have concentrated on explicit presentation of mathematical aspects of a particular application of 
non-Hermitian free random variables calculus to models which have the form of sums involving free unitary random 
matrices. We have derived general identities and solved three models as examples. 

The main aim of the paper is to convince the reader that the quaternion method offers a strikingly simple machinery 
to deal with certain instances of non-Hermitian models, reducing their solutions to elementary algebra. Sums of large 
free unitary matrices, or unitary and Hermitian, look at the first sight very complicated, but the quaternion technique 
provides a way to deal with them which involves only basic algebra. 

One question that can be posed deals with possible generalizations of the method to ensembles other than Hermitian 
and unitary used as "building blocks" of non-Hermitian models. This will depend on possible evaluation methods of 
B24J1 : in particular, this is doable when X satisfies (|26f) . 

One can also ask about some physical applications. The problem is that unitary matrices are not usually being 
added, but rather multiplied. The authors cannot see any applications at the moment. However, we see the solutions 
of sec. Mil 17x71 and Ivl as at least mathematically intresting, because we can precisely trace the way the unitary noise 
acts. For example, for the CUE plus GUE model, we would expect that the eigenvalues of the CUE matrix, lying on 
the unit circle, get smeared somehow when we impose the GUE noise; and this is exactly confirmed by our solution, 
which shows how the initial circle is smeared into the area between an internal circle and external ellipse. The three 
solutions we presented are thus quite intuitive, however we emphasize that the mathematical structure of the models 
is nevertheless involved, and this is the quaternion technique that allowed to tackle them. 

Another point may be to consider the general free additive unitary diffusion, which means the same infinite sum, 

E/r + ... + [% (12Q) 



M 

as in subsec. llV CI but for identically distributed arbitrary free unitary random matrices Ui, where eventually we are 
interested in the limit M — ► oo. This requires taking the equations (|74|l and (|75|l . together with general (|45|l . and 
plugging into them the expansions 

a = a ° + ~Trr ai + T7 fl2 + ■ ■ ■ ' b = b ° + ~7f7 bl + T7 &2 + ■ ' • ' ( 121 ) 
v M m V M 
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where we also write the standard Green's function for U as moment expansion, 

G u ( Z )= 1 - + ^ + ^ + .... (122) 

z z 2 - z 6 

Now comparison of appropriate terms in large- M expansions of the formuale will give some equations for ao, a\, . . . , 
and 6 0j .... Unfortunately, to our present knowledge, these equations are quite complicated and it may require a 
lot to extract the solution. This result will be an additive analog to the computation [jjfl of the free multiplicative 
unitary diffusion. 
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